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The lattice Boltzmann method with enhanced collisions and rest particles is used to calculate 
the flow in a two-dimensional lid-driven cavity. The abilitity of this method to compute the velocity 
and the pressure of an incompressible fluid in a geometry with Dirichlet and Neumann boundary 
conditions is verified by calculating a test-problem where the analytical solution is known. Different 
parameter configurations have been tested for Reynolds numbers from Re = 10 to Re = 2000. The 
vortex structure for a more generalized lid-driven cavity problem with a non-uniform top speed has 
been studied for various acpect ratios. 
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I. INTRODUCTION 

Lattice gas automata (LGA) are a rather new tech- 
nique in the world of fluid mechanics, but since the 
■ first discovery of a LGA by Frisch, Hasslacher and 
-»^Pomeau [1], which reproduces all terms of the Navier- 
^^Stokes equations, they have undertaken a fast develop- 
,-Hment. With the introduction of the lattice Boltzmann 
,__imethod (LBM) by McNamara and Zanetti [2] the sta- 
^tistical noise, present in the LGA, has been removed. 
'Higuera, Succi and Benzi developed the LBM with 'en- 
OOhanced collisions', where the foremost complicated col- 
lision operator was replaced by a much simpler one [3]. 
^-hA further simplification can be achieved if the collision 
^^erm is described by a simple relaxation. The LBM be- 
^^comes then a BGK model Qian [4]. The introduction of 
|; , - rest particles by H.Chen, S.Chen and Matthaeus [5] re- 
~l moved the unphysical factor in the pressure term of the 
Qvprevious models and allowed the correct calculation of 
-^ — .the pressure. 
^ In this paper the LBM with enhanced collisions and 
(3^ith rest particles is used to study the flow and pres- 
' sure distribution in a lid-driven cavity. There exists a 
Sfiarticular benchmark problem introduced by Shih, Tan, 
Gand Hwang [6], where the analytical solution is known. 
OThis gives the possibility to test the numerical code for 
^^oth Dirichlet and Neumann boundary conditions. Hav- 
ing verified that this method reproduces the analytic so- 
lution with small errors comparable to other numerical 
methods I go beyond the theoretical solvable problem and 
change the aspect ratio of the system. I study the change 
of the vortex structure during the increase of the aspect 
ratio for a small and a large Reynolds number [Re = 50 
and 1000). The occurring inflexional shear between the 
flows of opposite direction remains stable for Re = 50 
while for Re = 1000 it is unstable. 



II. THE LATTICE BOLTZMANN MODEL 

In general a lattice gas automaton consists of particles 
moving on links [i = 1, .., Mj) of a lattice from one node 
f%, to the next, where they can collide with one another. 
The time evolution for the population Ui for particles 
with speed cl is then given by a moving step and a colli- 
sion, described by the operator Qij:^ 

ni{ft: + Ci,i* + 1) = ni(T%,,i*) + Uij (1) 

Both operations, i.e. move and collision, must conserve 
mass and momentum. If the mean population are de- 
noted with Ni[fif,tif) =< ni[fif,tif) > the density p and 
velocity u is given by p = ^^ Ni and u = - ^^ NiCi, re- 
spectively. 

In the continuum limit the model should converge to- 
wards the Navier-Stokes equations. It has been shown 
that in two dimensions (2d) the hexagonal lattice is a 
proper choice [1], and in three dimensions (3d) it is the 
face-centered hypercubic (FCHC) lattice [7]. Assuming 
that all particles have the same speed and obey Fermi 
statistics'^ the equilibrium distribution for the population 
Ni is given by: 

N^ = [l + e''+^"' 



where h and g*are Langragian multipliers associated with 
the conserved quantities mass and momentum. For small 
Mach numbers one can compute h and q pertubatively 
in u and end up with the following expansion for .AT?': 

7V/'= d+d—{ci)cUc+d — g{p)QiapUo,up +0(u^) 
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and Qiap = {ci)a{ci)p - J^^ap 

d is the density per link d = -^. 

The continuum limit is performed as a multiscale expan- 
sion in a small parameter e, which can be identified with 
the local Knudsen number. The actual population Ni is 
expanded around its equilibrium: 

Ni = N!'' + eNr"'^ + Oie'') 



The time step is set equal to 1. 

The Fermi statistic has its origin in the lattice gas mod- 
els with Boolean variables. Within the lattice Boltzmann ap- 
proach it is possible to choose a different distribution function, 
see e.g. [4]. 



Using the mass and momentum conservation for the dy- 
namics of the lattice automaton one arrives at the Navier- 
Stokes equations [1]: 



dp dup 
dt drp 



(4a) 
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+ 0{eu^) (46) 

t = etif + e^tif and f = ef%, are the continous time and 
space varibiable, respectively. The viscous stress tensor 
Sap is given by: 

fd{pup) d{pua) 2 d{pu^ 

^nfl = V —On ' 



\ dra drp D dr-y 

v is the viscosity of the system, which depends on the 
details of the collisions. An expression for v will be given 
below. The momentum flux tensor H-ap has the following 
form: 

^<xp = pc](l-9{p)-^\+P9{p)uaUp (6) 



Cj = ^{c? jT)) is the sound speed and p the pressure. 
One disadvantage of the lattice automaton with Boolean 
variables is the need of performing a statistical average 
numerically. This step can be made obsolete by using the 
mean population TV^, which defines the lattice Boltzmann 
model (LBM). If the Boltzmann approximation is taken 
for the collision, the collision operator 0^^ can be approx- 
imated by an expression linearized in the non-equilibrium 
part 7\r"°' = Ni — 7\r?' of the population numbers: 

^ij ^ Aij{Nj - N;^) (7) 

The elements aij of the matrix Aij describe the scatter- 
ing between direction i and j and depend on the angle 
between the directions only. This approach was intro- 
duced by Higuera, Succi and Benzi and is called LBM 
with enhanced collisions. Because Aij is cyclic and sym- 
metric, their elements aij can be expressed in terms of 
the three non-zero eigenvalues of the matrix: A, a and 
T [3]. The influence of a and r on the results is negligi- 
ble, and they both can be set equal to -1 [3]. The other 
eigenvalue is related to the viscosity via: 




Furthermore, it can be shown that — t is the relaxation 
time, in which the population Ni at a node converges to 
its equilibrium N^'^ [8]. 

For the one-speed model, described above, the pressure 
p depends explicitely on the velocity u, as it can be seen 



from equation (6). This unsatisfying feature can be re- 
moved by introducing a reservoir of rest particles with a 
density do which is a function of the density d [5]. The 
total density of the system is then given by: 

p = do + Msd (9) 

In the following I will restrict myself to the FCHC lat- 
tice where D = 4, Mj = 24 and c = 2. In this model 
choose the density of the rest particles is chosen as: 

Formally, this fixes the ^-factor to 1. But because the 
population of rest particles is local and not moving, it can 
be chosen any other density do with the only restriction 
that this density is large enough to activate the rest par- 
ticles for the highest velocity. Since it will be seen later 
that also the density d has no influence on our results, 
both d and do can act as parameters for additional phys- 
ical modelling, e.g. in multi-phase or multi-component 
flows. 
The equilibrium distributions can now be written as 

TV/' = d+ — 1^2(0^)^^^ + 3{ci)c{ci)pUcUp - u^J (11a) 

for the moving particles and 
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^' = '^0-^ (116) 

for the rest particles. The population of the particles No 
relaxes into its equilibrium with the relaxation time of 
the system, i.e. with T^ei = —j- Therefore, the time 
evolution of the population of rest particles is given by: 

No{f,,U + 1) = No{f,,U) + \{No - <') (12a) 

Because the rest particles are activated isotropically into 
moving particles, the time evolution of the populations 
of the moving particles has the following form: 

Ni{f, + Ci,U + 1) = Ni{f,,U) + Aij{Nj - 7V/«) 

-^{No-N^') (126) 

In this paper I use a projection of the FCHC lattice on 
two dimensions which is described in detail in [8]. In total 
there are nine different populations for the moving parti- 
cles, four in the axial directions [i = 1, 3, 5, 7) with Ci = 

i cos ^-^-TT, sin ^-^-TT ) which have a weight of 4 due to the 

fact that each is the projection of four links from the orig- 
inal space, four in the diagonal direction [i = 2, 4, 6, 8) 

with Ci = ( cos ^^TT, sin ^^TT J , and one perpendicular 

to the plane, which has also a weight of 4. The tenth 
population is that of rest particles. Therefore, the total 



density is given by p = ^ TV^ + 4 ^ TV^. Introducing 

i even i odd 

new populations via Nl=Ni — d and Nq=No — do and 
remembering that, due to the symmetry of the collision 
matrix, Aij^Ni—N^'^) can be replaced by Aij^Ni — N^'^' ) 
[8], equation (10a) and (106) can be rewritten in a form, 
which is more convenient for computational purposes: 

N'(K,U) = N'(K,U -1) + \{N' + ^] (13a) 
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N!{f^,U) = N!{f^ - Ci,U - 1) + Aij 



N-{f^,t*-l) - -{ci)a{ci)pUcUp + —u^ 
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Instead of calculating the total density we now com- 
pute bp= Yii ^i +^ Yl ^i- The velocity is still given 

i even i odd 

by expression u = - ^^ Ni Ci . 

The hydrodynamical system is described by the 
Reynolds number Re = -'^^^ with the characteristic ve- 
locity uq and the number of mesh points Iq for the char- 
acteristic length. As a consequence, two parameters out 
of i^, Iq, or uq can be chosen while the third is fixed by 
the Reynolds number. The used LBM implies some re- 
strictions to the velocity, where an upper limit of u = 0.2 
(Mach number M = 0.28) has been found [9], and to 
the viscosity where lower limit of i^ = 0.14 has been ob- 
served [9]. The upper limit for velocity results from the 
compressibility of the used lattice gas, because the spatial 
change of the density describes the pressure (see equation 
(6)). If the velocity u becomes too large the compressibil- 
ity effect is not negligible any more. The lower limit for 
the viscosity has its origin in the relaxation behaviour of 
the system, because —A = g-^rr is the relaxation param- 
eter uj and tends to 2 if the viscosity is decreased towards 
zero. If the overrelaxtion (w > 1) is too large the system 
becomes unstable. The lower limit for ly mentioned above 
was derived from the numerical experience with the LBM 
without rest particles. The behaviour of the LBM with 
rest particles is discussed in section III. 

From now on I change the notation in the following 
way: parameters in the LBM are labeled with lb to dis- 
tinguish them from the values in the 'real world'. I define 
u and I via 



U = Ui^B ■ u 



and I = ZiB • I (14a), 



respectively. If the normalized pressure p = ^ is used, 
where p' is pressure and p the density, the normalized 
pressure can be obtained from the variation ^p^B in the 
mean density p^e via: 

P=- — (146) 

Plb 2 

The relation for a body force / is: 

/ = y.B (14c) 



In the following I will give a short description how I 
implement the boundary conditions. I impose the bound- 
aries on half way between two nodes. At the left, right, 
and bottom boundary there are no-slip conditions and 
particles are bounced back. This takes place for all Di- 
rechlet boundary conditions except that for non-zero ve- 
locities one has to add the right momentum to the par- 
ticles travelling on the diagonal links. Considering the 
top boundary, the population N3 is copied into Nr with- 
out adding an extr a m omentum. The population N2 is 
copied into Ne and ^^u/b is subtracted, u/b is the (lo- 
cal) velocity at the top given by the boundary condition. 
A similar process occu rs to the population N4, which is 
converted into Ns and ^^Uj°b is added. 
In the case of Neumann boundary conditions particles 
travelling on diagonal links undergo a specular reflec- 
tion. The required derivative at the boundary 5uj°b is 
achieved by adding (or subtracting) the moment loss due 

to viscous stress ^i^^B<5uI°| = ^i^^B^^^f^f^ 



to 

top 



the reflected particles. 
At this point I have to make a comment about the den- 
sity used in equation (13a) and (13b) and in the update 
at the boundary. Within the numerical errors I do not 
observe a difference between using the local density p^B 
and using the mean density p^B in (13a) and (13b) for 
the calculations of the test-problem (see next paragraph). 
This justifies to use the mean value also for the update 
of the boundary conditions. Using p^B instead of the lo- 
cal PiB for calculating the equilibrium distribution saves 
computer time in the collision step. 

The great advantage of the LBM's is the local nature of 
the boundary conditions. Though it can be shown nu- 
merically that the Neumann boundary condition imple- 
mented in the described way is of second order in space, 
only the nearest neighbours are involved in its calcula- 
tion. 

After all the theoretical considerations one ends up 
with a very simple algorithm consisting out of three sub- 
steps per time-step: 

1. The boundary conditions are set according to the 
previous time step. The populations N^ are calcu- 
lated at the boundary. 

This step involves nearest neighbours. 

2. At every node new populations are calculated due 
to the moving of its nearest neighbours. 

This step involves nearest neighbours. 

3. At every node the collision term is calculated and 
the body forces are added. 

This step is completely local. 

For all calculations I use the zero-speed initial condition, 
i.e. N;, = N! = 0. 

I run the code on a DEC Alpha station for lattices up 
to 150 X 150 meshpoints and on a NEC SX-3-24R for 



larger lattices. On the DEC about 10/xs are needed per 
timestep and node, on the NEC about 83ns. On the NEC 
I achieve a performance of between 2.0 and 2.4 GFlops, 
which has to be related to a peak performance of 6.4 
GFlops. The program has not been specially adapted for 
the NEC computer. 



III. TEST-PROBLEM 

To test the LBE method with the previous described 
boundary conditions I choose a benchmark problem orig- 
inally proposed by Shih, Tan, and Hwang [6]. The geom- 
etry of the problem is a square box of Imx Ira [xi^n = I'n 
and yi^n = I'n) with no-slip boundary conditions at the 
left, right, and bottom boundary and a shear flow at the 
top with the velocity profile: 

■u^(a;) = 16(a;^-2a;^ + a;2) (15) 

The maximum speed of Ira/ s is chosen as the characteris- 
tic velocity and therefore the Reynolds number becomes: 

1 m^ 

Re = (16) 

V s 

Assuming a time-independent solution of the Navier- 
Stokes equations, which is true for small Reynolds num- 
bers [Re < 1000), the steady state velocity u°° = 



°° U°°] 
■x J "-y ) 


is set to: 










u^= 8s{x)v'{y) 


(17a) 






u^ = -8s'{x)v{y) 


(176) 




with 


s{x) = E^ - 2x^ + x^ 






and 


■"{y) = y^ -y^ 





If the pressure is given by: 

p°° = ^ • 8{S{x)v"'{y) + s'{x)v'{y)) + 6i{S2{x){v{y)v" {y) - [v'iy))') (18) 

the solution of the time-independent Navier-Stokes equations results in a body force / = [f^, fy): 

U = (19a) 

fy=iy 8{24:S{x) + 2s'{x)v"{y) + s"'{x)v{y)) + 64:{S2{x)V{y) - v{y)v' {y)Si{x)) (196) 

The above used functions are defined as follows: 

S[x) = J s[x)dx 

Si{x) = s{x)s"{x) - {s'{x)Y 

S2{x) = J s(x)s' (x)dx 

V{y) = v{yy"{y)-v'{y)v"{y) 

I start with a detailed study at Re = 50 to check the 
influence of the density p^B, the (characteristic) velocity 
UiB and the lattice size on the accuracy of the results. 



The difference between the theoretical and the numer- 
ical results is characterized by the relative error of the 
velocity and the pressure 
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and 



TEgHz? (m) 

respectively. I do not find an influence of the density p^e 
or the relation do/d on the results. The errors for both 
velocity and pressure become smaller with increasing lat- 
tice size and decreasing velocity as it is shown in Fig.l. 
For the fixed velocity u^b = 0.02 the error decreases 
approximately linear with the mesh size, for higher ve- 
locities the slope is even less. A linear behaviour was 
also observed for calculations with the LBGK method 
{ulb = 0.1) [10]. For a fixed mesh size the axis of veloc- 
ity can be reinterpretated as the timestep Ai in the 'real 
world', because the timestep is given by Ai = ULs/ly 
where ly is the number of mesh points in the vertical. 
The required computer time increases linearly with de- 
creasing velocity U£,b, but it is proportional to the third 
power of the lattice size. Therefore, one would prefer to 
reduce the velocity rather than increasing the lattice size 
but here one is limited to (ui^B)rnin = °';''™'" . From 
Fig.l one can deduce that there exist a minimum for the 
error of the velocity at a certain {ui^B)min- Numerically 
it is found that this value is smaller the larger the lattice 
is. The behaviour of the pressure is slightly different und 
the minimum lies at higher values of u^b. Going beyond 
this value, irregularities in the pressure field grow in the 
left upper corner. As an example the results of a calcu- 
lation for i^iB = 2.5 • 10~^ and ly = 50 for Re = 100 are 
shown in Fig. 4a. Though the error of the pressure field 
is rather small (see Table I), some abnormal behaviour 
of the pressure is observed. The velocity field instead 
exhibits no unusual featurs compared to those of Fig. 3 
and it is still very close to the analytical solution. Using 
the original LBM without rest particles the origin of the 
vortex is shifted (see Fig. 4b) and the calculation loose its 
validity.^ For larger values of the viscosity I do not find 
a difference between the two kinds of methods within 
the numerical errors. This means that the reservoir of 
rest particles stabilizes the system for short relaxation 
times. With the LBGK method including rest particles 
Hou et.al. reach a lower limit of the viscosity of about 
i^LB = 2.56 • 10~^ [10]. For their calculations with a fixed 
top velocity of u^b = 0.1 they used a 256 x 256 lattice for 



In this model the pressure cannot be calculated directly via 
equation (14b). 



all Reynolds numbers so that for e.g. Re = 1000 the vis- 
cosity is i^ = 0.0256. There is no possibility to compare 
their results for the pressure with other numerical calcu- 
lations nor did they publish a comparison of the pressure 
field with different values for the viscosity by changing 
the lattice size. 

For the Reynolds number Re = 100 I perform three 
calculations with lattice sizes 50 x 50, 100 x 100, and 
150 X 150. The behaviour of the errors in time is shown 
in Fig. 2a. Since the momentum is mainly transported by 
diffusion, the error of the velocity decreases exponentially 
in time and the exponent is proportional to the viscos- 
ity. Therefore, one expects a similar time-behaviour of 
the error for Re = 10 on a time scale of one magnitude 
smaller. Indeed, a calculation for Re=10 with a lattice 
size of 50 X 50 and Ui,B=0.04 shows this behaviour of the 
error (see Fig. 2b). The oscillatory behaviour of the er- 
rors can be reduced by decreasing both the velocity u^b 
and the viscosity i^^e keeping the lattice size fixed. The 
reason for the oscillitory behaviour might be the long 
relaxation time (high viscosity) and the large timestep 
(high velocity), which may result in a problem to relax 
the system into its local equilibrium in every timestep. 
In Fig. 2b the errors of my calculations are compared with 
those of a finite element calculation with a fractional step 
scheme where the lattice size was 64 x 64 [11].'* 

I go beyond the Reynolds number of 100 and perform 
calculations for Re = 500, 1000, 2000, and 5000. For 
Re = 5000 computation becomes extremely time con- 
suming because of the large lattice size. I try to reduce 
the lattice size to 1000 x 1000 by choosing three different 
parameter configurations: a) u^b = 0.1 and i^^e = 0.02, 
b) UiB = 0.05 and i^^e = 0.01, and c) u^b = 0.025 
and i^iB = 0.005. In the first two cases a blow up of 
the density in the upper right corner takes place after 
some thousand iterations. While in these both cases the 
timestep in real units is larger than for the calculation for 
the 2000x2000 lattice (a) Ai = MO'^, b) Ai = 5-10"^) 
it is the same in the third case (Ai = 2.5 • 10~^s). In the 
third case due to the small viscosity unphysical pressure 
fluctuations arise on the right-hand side and lead to prob- 
lems calculating the collision term. From the calculations 
with the 2000 x 2000 mesh I observe a deminishing of the 
vorticity with time but it is not yet clear if the system 
will converge towards a steady-state like solution with 
only small perturbations in time, or if the center of the 
main vortex will oscillate and some extra vorticities will 
remain. I have to stress that this case is beyond the va- 
lidity of the analytical solution and the behaviour of the 
system is not yet known. 

Instead of applying the Dirichlet boundary condition 
at the top one can force the flow by the Neumann bound- 



'Only the error for the velocity was available. 



ary condition: 

du^{x,y) 



dy 



SOfx"^ - 2x^ + x^) (21) 



I test this boundary condition for Reynolds numbers 
Re = 50, Re = 100 and Re = 1000. The results are 
listed in Table II, which shows that the Dirichlet and the 
Neumann boundary condition give results of nearly the 
same accuracy. The convergence rate for the errors of 
the calculation for Re = 100 can be seen in Fig. 2a. It 
is of the same order as for the Dirichlet boundary condi- 
tion. The fact that the LBM yield the result regardless 
of the kind of boundary condition is very striking in the 
sense that the numerical effort for both boundary con- 
ditions is the same while for the classical methods it is 
larger for the Neumann boundary condition than for the 
Dirichlet. The very good agreement of our numerical re- 
sults with the analytical solution encourage me to use 
the Neumann boundary condition later for the extended 
cavity problem. 



IV. EXTENDED DRIVEN CAVITY PROBLEM 

In the previous section I have verified that the LBM 
with enhanced collisions and rest particles reproduces the 
flow in a lid-driven cavity very nicely. Now I turn to study 
the flow in a cavity without applying the body force and 
extend our studies to aspect ratios Ar := ^^^^ different 
from 1. I still apply the same Dirichlet, or in a second 
step the same Neumann boundary condition at the top 
of the cavity as in the test-problem. The two different 
boundary conditions have no longer the same meaning 
because velocity field is now unknown and not given by 
equation (17). I have to be a little bit more specific about 
the geometry of the system and the boundary conditions: 
the y-axis is fixed to length 1 for all calculations [yien = 
Im). Defining a new variable x = Ar~^ ■ x, with < £ < 
1, the Dirichlet boundary condition is written as 

■u^(a;, 1) = 16(£^-2£^ + £2) (22) 

and the Neumann boundary condition as 



du^{x,y) 



dy 



80(£^-2£^ + £2) (23). 



Therefore, in the case of the Dirichlet boundary condition 
the maximum velocity at the top is still Ira/ s (and the 
average velocity 0.533?n/s) but the velocity gradient in 
the x-direction decreases as the aspect ratio increases: 

du^{x^ 1) ^ 3g^^-i(g-3 _ 3-2 ^ -) (24) 

In other words the velocity field becomes smoother for 
larger aspect ratios. The same is true for the Neumann 
boundary condition. 
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I perform runs for two different values of the vis- 
cosity {ly = 0.02m?s~^ and ly = 0.001m~^s~^, which 
correspond for the Dirichlet boundary condition to the 
Reynolds numbers Re = 50 and Re = 1000) and for four 
different aspect ratios [Ar=0.2, 1, 10 and 50). I also 
perform calculations with the Dirichlet boundary condi- 
tion at those Reynolds numbers, which I observe with 
the Neumann boundary condition and ly = 0.001?n~'^s~-'^ 
at the particular aspect ratio. Table III lists up all the 
calculations with the used parameters and the origin of 
the vortices. 

I start discussing the results for the small Reynolds 
number. Without the body force the velocity field is no 
longer symmetric according to the plane at a; = 0.5?n 
but the origin of the vortex is now at f = (0.56, 0.76)?n 
instead oi f = (0.5, 0.707)?n for the test-problem. With 
increasing aspect ratio the vortex is stretched and it is 
not very helpful to define an origin of a vortex in this 
case. For Ar = 50 the velocity profile for u^ is parabolic 
except in the vicinity of the left and right wall. The in- 
flection line between forward and back flow lies very near 
the theoretical value [y = 2/3yie„) for an infinitily long 
cavity with a constant flow (constant top velocity). In 
Fig. 5 the velocity profile of u^ is shown at vertical cuts 
through the cavity. The deviation of the velocity at the 
top and the resulting fluid exchange in the vertical is too 
small to disturb this flow structure. 

A different behaviour of the flow with increasing aspect 
ratio is observed for a Reynolds number of Re = 1000. 
For the aspect ratio Ar = 1 the origin of the vor- 
ticity is observed near the centre of the cavity [f = 
(0.54, 0.573)?n) as it is shown in Fig. 6. The main driv- 
ing force of the system is now the convective momentum 
transfer. The deep pressure at the centre balances the 
Centrifugal force of the strong circulation. Pressure and 
velocity field are quite similar to those of the driven cav- 
ity problem with a fixed top velocity [10]. 

The increase of the aspect ratio to Ar = 10 does not 
lead to just a stretching of the vortex but a vortex is 
created near the right-hand side of the cavity (its origin 
is at F = (8.84, 0.522)?n ). The vorticity of this vortex 
is much stronger than those of the vortex present in the 
case of Ar = 1, because the main part of the energy, 
which is put into system by the forced flow at the top, 
goes into this vortex. The energy of the vortex is so high 
that the fluid leaving the vortex has enough kinetic en- 
ergy to run against the pressure gradient of the flow (see 
Fig. 7). The situation changes totally if the aspect ratio 
is further increased to a value of Ar = 50. The vortex 
on the right-hand side has now disappeared (see Fig. 8). 
I presume that the deviation of the velocity at the top 
is too smooth to cause a vortex structure. Nevertheless, 
the system is unstable in the sense that pressure waves 
are observed travelling in the horizontal direction. The 
origin of these waves is not yet clear and is a topic of 
recent research. Another part of the future studies is the 
question, at which aspect ratio does the transition from 
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a structure with a vortex to a vortex-free structure oc- 
cur. The results of these investigations will be published 
elsewhere. 

I now turn to present the results for the Neumann 
boundary condition at the top and the low viscosity. 
For an aspect ratio of Ar = 1 a maximum top speed 
of Umax = 0.295?n/s is found, which means that the 
Reynolds number is Re = 295. Computation of the pres- 
sure and velocity field with the Dirichlet boundary condi- 
tion at this Reynolds number shows that the type of the 
boundary condition has nearly no influence on the result- 
ing pressure and velocity (see also Table III). The small 
shift of the maximum top-speed towards the right-hand 
side in the case of the Neumann boundary condition has 
no effect on the structure of the vortex. The situation 
changes if the aspect ratio is Ar = 10. In Fig. 10 two vor- 
tices are found on the right-hand side, a large one similar 
to those for the Dirichlet boundary condition with its ori- 
gin at f = (9.16, 0.54)?nand a small one at the upper left 
part of the large one with its origin at f = (7.59, 0.63)?n. 
In accordence with the difference in the velocity field the 
pressure field exhibits now two high pressure cusps in- 
stead one for the Dirichlet boundary condition. In this 
case the maximum speed is Umax = 0.58m/ s correspond- 
ing to a Reynolds number of Re = 580. For the calcula- 
tions with the Dirichlet boundary condition no difference 
in the flow and pressure field on principal can be found 
between Re = 1000 and Re = 580 (Fig. 11). In both 
cases a high pressure zone is build up at the right-hand 
side excluding the main flux from this area. This area is 
missing for the Neumann boundary condition. 
One possibilty to increase the Reynolds number is in- 
creasing the deviation of u^ at the top, given by equation 
(24), by a certain factor. If a factor 2 is chosen a Reynolds 
number Re = 830 is observed with a flow and pressure 
pattern shown in Fig. 12. The second high pressure cusp 
has disappeared nearly completely and the smaller vortex 
is shifted towards the inflection line between the flows of 
opposite direction. In Fig. 13 the dependency of the top- 
velocity on the horizontal position is presented for differ- 
ent Reynolds numbers. For high Reynolds numbers the 
dependency is quite different from that for the Dirichlet 
boundary condition while for the low Reynolds number 
only the center is shifted. 

If the aspect ratio is increased to Ar = 50, the velocity 
field is qualitively the same as for the Dirichlet boundary 
condition. The Reynolds number is Re = 880. This case 
exhibits also pressure fluctuations. 

Neumann boundary conditions are of great interest in 
the case where the fluid in the cavity is coupled to a mov- 
ing gas at its top. For example, in shallow water simu- 
lations the derivative at the top is proportional to the 
velocity of the wind. In Fig. 13 the Neumann boundary 
condition on top of the box is plotted with the velocity 
and pressure field, and this distribution can be regarded 
as the velocity of the wind up to a certain factor. In the 
future the calculations will be extended to three dimen- 
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sions. In principal, there should be no difference in the 
numerical behaviour with respect to the two-dimensional 
case because the method, used in this paper, is a projec- 
tion of the three-dimensional version. 



V. CONCLUSION 

In this paper I have shown that the lattice Boltzmann 
method with enhanced collisions and rest particles is a 
fast and accurate method for calculating the velocity and 
the pressure in a driven cavity by comparing our numer- 
ical results with the analytical solution of the problem. 
Using the Dirichlet or the Neumann boundary condition 
for the forced flow at the top makes no difference ac- 
cording to the convergence of the numerical results. This 
is a very striking observation because the implementa- 
tion of the Neumann boundary condition in the LBM 
is extremely easy and involves only nearest neighbours. 
It is possible to decrease the viscosity beyond the value 
v^^ = 0.014, which was found to be the lower limit for 
the LBM without rest particles. The pressure is more 
sensitive to relaxation problems than the velocity. 
In the second part of this paper I have presented the 
results for a driven cavity problem with a non-uniform 
top speed and aspect ratios different from one. For small 
Reynolds numbers [Re = 50) the result of increasing the 
aspect ratio is a stretching of the original vortex. At 
higher Reynolds numbers [Re = 580, 1000) a strong vor- 
tex on the right-hand side of the cavity has been observed 
for an aspect ratio of ^t" = 10. At this aspect ratio a sig- 
nificant difference has been found between the Dirichlet 
and the Neumann boundary condition because a second 
but smaller vortex does exist if the Neumann boundary 
condition are applied at the top. 
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TABLE I. Parameters and errors for the test-problem with 
the Dirichlet boundary condition at the top. At is the 
timestep of the calculation in terms of the 'real world'. The 
total time indicates the time in which the system has been 
converged into its steady-state within the numerical errors 
except for Re = 5000 where the calculation is extremely time 
consuming. 

The CPU-time shows the consumed computer time on a 
DEC- Alpha station (a) or a NEC SX-3-42R (SX-3). 
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TABLE IL Parameters and errors for the test-problem 
with the Neumann boundary condition at the top. At is the 
timestep of the calculation in terms of the 'real world'. The 
total time indicates the time in which the system has been 
converged into its steady-state within the numerical errors. 
The CPU-time shows the consumed computer time on a 
DEC- Alpha station (a) or a NEC SX-3-42R (SX-3). 
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TABLE III. Parameters and results for the extended driven 
cavity problem. 

Re is the Reynolds number, Ar the aspect ratio, be denotes 
the boundary condition (D: Dirichlet, N: Neumann), f are 
the coordinates of the origin of a vortex, and Umax is the 
maximum speed (allways at y = Im). 
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FIG. 1. Relative errors for velocity and pressure according 
to the definition of equation (20) for Re = 50 and different 
parameter configurations. The axis for the velocity Mi,b is 
logarithmic for a better visualisation of the results. The lines 
on the basic plane mark contours of constant viscosity Ulb- 
The values vor Ui^b are from left to right: 0.0025, 0.005, 0.01, 
0.05, 0.1, and 0.2. 



FIG. 2. Relative errors for velocity and pressure according 
to the definition of equation (20) versus the time in the 'real 
world' for Re = 100 (left) and Re = 10 (right). Lines with 
points in its style belong to the pressure, the others to the 
velocity. 

For Re = 100 all results as listed in Table I are plotted except 
the case with the viscosity u = 0.0025. Also the result with 
the Neumann boundary condition is plotted (compare Table 

11). 

For Re = 10 the results for Mi,b = 0.04 (solid line for velocity 
and pointed line for pressure) and Mi,b = 0.004 (dashed line 
for velocity and dashed-pointed for pressure) are shown. 
The line with the stars as markers is the result of a finite 
element calculation [11]. 



FIG. 3. Pressure and velocity profile for the test-problem 
{Re = 100). 



FIG. 4. Pressure and velocity profile for the test-problem 
[Re = 100). The viscosity is vi,b = 0.0025. 



FIG. 5. Velocity profiles for Ux at different vertical cuts 
through the cavity for Re = 50 and Ar = 50. 



FIG. 6. Pressure and velocity profile for Re = 1000, aspect 
ratio Ar = 1 and the Dirichlet boundary condition at the top. 



FIG. 7. Pressure and velocity profile for Re = 1000, aspect 
ratio Ar = 10, and the Dirichlet boundary condition at the 
top. 

The aspect ratio of the picture does not agree with those of 
the calculation. 



FIG. 8. Velocity profile for Re = 1000, aspect ratio 
Ar = 50, and the Dirichlet boundary condition at the top. 
The aspect ratio of the picture does not agree with those of 
the calculation. 



FIG. 9. Pressure and velocity profile for Re = 295, aspect 
ratio Ar = 1, and the Neumann boundary condition at the 
top. 
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FIG. 10. Pressure and velocity profile for Re = 580, aspect 
ratio Ar = 10, and the Neumann boundary condition at the 
top. The aspect ratio of the picture does not agree with those 
of the calculation. 



FIG. 11. Pressure and velocity profile for Re = 580, aspect 
ratio Ar = 10, and the Dirichlet boundary condition at the 
top. The maximum speed at the top is reduced to 0.58m/s, 
the viscosity is j/ = O.OOlm /s. 

The aspect ratio of the picture does not agree with those of 
the calculation. 



FIG. 12. Pressure and velocity profile for Re = 830, aspect 
ratio Ar = 10, and the Neumann boundary condition at the 
top. The maximum derivative at the top is increased to 10s~ , 
the viscosity is j/ = O.OOlm /s. 

The aspect ratio of the picture does not agree with those of 
the calculation. 



FIG. 13. Velocity at the top for different Reynolds numbers 
when the Neumann boundary condition is applied at the top. 
The aspect ratio is Ar = 10, the viscosity v = O.OOlm /s for 
the high Reynolds numbers, and v = 0.02Tn /s for the low 
Reynolds number. 
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